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ABSTRACT 

The effect of gas ejection on the structure and binding energy of newly formed steUar 
clusters is investigated. The star formation efficiency (SFE), necessary for forming a 
gravitationally bound stellar cluster, is determined. 

Two sets of numerical N-body simulations are presented: As a first simplified 
approach we treat the residual gas as an external potential. The gas expulsion is ap- 
proximated by reducing the gas mass to zero on a given timescale, which is treated 
as a free parameter. In a second set of simulations we use smoothed particle hydrody- 
namics (SPH) to follow the dynamics of the outflowing residual gas self-consistently. 
We investigate cases where gas outflow is induced by an outwards propagating shock 
front and where the whole gas cloud is heated homogeneously, leading to ejection. 

If the stars are in virial equilibrium with the gaseous environment initially, bound 
clusters only form in regions where the local SFE is larger than 50% or where the 
gas expulsion timescale is long compared to the dynamical timescale. A small initial 
velocity dispersion of the stars leads to a compaction of the cluster during the expulsion 
phase and reduces the SFE needed to form bound clusters to less than 10%. 

Key words: globular clusters: general - open clusters: general - stars: formation. 
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1 INTRODUCTION 

During the formation of star clusters gas expulsion, caused 
by feedback of young massive stars, terminates the star 
formation epoch and can unbind the stellar system. Many 
mechanisms leading to gas loss exist, like ionizing radiation, 
stellar winds or supernova explosions. It is still uncertain 
which of those play the major role. 

The gas expul sion reduces the binding energy of the 
cluster. Hills (19801) showed, using analytic approximations, 
that a system of stars and gas loosing more than half of its 
mass in less than a crossing time will disrupt. Thus, to ob- 
tain bound stellar clusters, a star formation efficiency (SFE) 
e = Ms/Mc > 0.5 is needed, where Ms and Mc are the mass 
of the stellar component and of the initial gas cloud prior 
to star formation, respectively. As the typical SFEs are less 
than 10%, the formation of gravitationally bound old open 
clusters and globular clusters is an interesting and yet un- 
solved problem. 

Numerical simulations investigating the stability of 
young star clusters after gas expu lsion have been done by 
Lada, Margulis & Dearborn (1984). They showed that open 
star clusters, initially in virial equilibrium with the sur- 
rounding residual gas and containing up to 100 stars, can 
remain bound even if the SFE is as low as 30%. In their 
simulations they treated the residual gas as a variable ex- 



ternal potential added to that of the stars. Goodwin (1997) 
extended these simulations to globular clusters (GC), in- 
creasing the number of particles, allowing for different gas 
expulsion mechanisms and including loss of st ars due to 
the galactic tidal field. Recently, Adams (200C) presented 
a semi-analytic model for the formation of bound star clus- 
ters even for global SFEs much smaller than 50%. 

We present a new set of numerical simulations of the 
early evolution of young GCs, covering a broad range of 
SFEs and gas expulsion timescales. Because the typical col- 
lision dominated relaxation timescale of GCs is much larger 
than the dynamical timescale, we restrict ourselves to colli- 
sionless N-body calculations for the stellar component. We 
present two sets of simulations: As a first assumption, in 
Sect, ti we treat the residual gas as an external potential. 
We investigate the dynamics of the cluster during and after 
gas expulsion and derive constraints for the SFE required 
to form a bound cluster. Using a combined N-body and hy- 
drodynamic code, in Sect, a we investigate the gas removal 
more self-consistently. Conclusions follow in Sect. H. 
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Table 1. Parameters of the initial configurations (N-body) 



model 


Wo 


tc 


Rt 


N 


5 


Nl 


3.0 


0.28 


1.26 


1000 


0.1 


N2 


5.1 


0.17 


1.33 


1000 


0.1 


N3 


3.0 


0.28 


1.26 


4000 


0.05 



Wo = ?/'(0)/(T^: scal ed cen tral potential of the King profile, see 
Binney & Tremaine (1987); tc'. crossing time at half-mass radius; 
Rf. tidal radius of King profile; A'^: number of particles used in 
simulation; S: numerical (Plummer) smoothing length; total mass 
(stars and gas) of all models was set to 1; all quantities are given 
in dimensionless code units (see text) 



2 GAS EXPULSION IN PURE N 
SIMULATIONS 



BODY 



Our simulations start after the cluster has formed, but be- 
fore the residual gas has been expelled. The stars are rep- 
resented by collisionless, equal mass N-body particles. The 
effect of the ejection of the residual gas on the stellar sys- 
tem is treated as a time variabl e ext ernal potential, similar 
to the approach of Lada et al. (1984). 



In the following, all quantities are given in dimensionless 
code units (gravitational constant G = 1). Thus, taking typ- 
ical globular cluster properties M = 10^ M© and R = 10 pc 
as mass and length units, respectively, we obtain a time unit 
i = [G p)~^'^ = 1.5 10® yr (equal to the dynamical or cross- 
ing timescale) with the density unit p = M / B? . 

The N-b odv calculations are done using a hierarchical 
tree method (Barnes fc Hut 1986). 



2.1 Initial Configuration and the Gas Expulsion 
Phase 

To obtain a stable initial configuration, in a firs t step the 
stars are distributed according to a King (1966) distribu- 



tion function with total mass equal to that of the initial gas 
cloud. The potential is tabulated and is used for modelling 
the external gas potential during the simulation. Therefore, 
stars and gas have equal density distributions. Finally, the 
mass of the stars and the gas are scaled according to the 
SFE. Now the stars are in virial equilibrium within the sum 
of their own potential and the potential of the gaseous com- 
ponent. The parameters of the different models are given in 
Table pi. For each model, the SFE is varied between 0.15 and 
0.80. 

To test whether the initial system is in virial equilib- 
rium, several calculations without gas expulsion are per- 
formed. The density distributions are well conserved. 

We simulate the gas expulsion starting at time t — to 
by multiplying the external potential by a time dependent 
factor 



(t-to)/to 



t < to 

to < t < to + tcxp 

t > to + toxp 



(1) 



where ioxp is the time that is needed to drive the gas out of 
the system (gas expulsion time). 

We can estimate the order of the gas expulsion time: 
The isothermal sound speed of a molecular cloud gas with 
temperature T = lOK and molecular weight fi — 2.36 is 
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. t=494 
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Figure 1. Time evolution of model N3. The time t is given in 
dimensionless units. The SFE is e = 0.4, expulsion time icxp = 2. 
The plots show the N— body particles (small dots) projected onto 
the x—y— plane; unbound particles are indicated by squares. 



a — y^Rgas T/fi = 0.19 km s^^, where Rg^n is the gas con- 
stant. If we consider a disruptive process that starts at the 
centre of a cloud as given by model Nl and travels out- 
wards with sound speed, it will need approximately a time 
of icxp ~ Rt o,~^ ~ 6.6 10® yr (or texp ~ 4, given the dimen- 
sionless code units above) to reach the edge of the cloud. Fast 
processes (e.g. supernova explosions) may remove the gas on 
shorter timescales. The gas expulsion time will therefore pre- 
sumably be of the order of a dynamical time which is equal 
to the unit of time. In the simulations we use texp = 0, 2, 4 
and 10, which are equal to 0, 7, 14 and 36 crossing times at 
half-mass radius of model Nl and N3 and 0, 12, 24 and 59 
crossing times of model N2. 



2.2 Dynamics of the Cluster During and After 
Gas Expulsion 

The typical evolution of the N-body part of a cluster with 
tcxp = 2 (7 crossing timescales) is displayed in Fig. hi. Start- 
ing at t = 20 the external potential is slowly reduced to zero 
as described in the previous section. The cluster expands. A 
certain amount of stars gets unbound and starts leaving the 
system. The bound ones relax after the gas has been com- 
pletely removed, forming a broader configuration. A particle 
is believed to be unbound if its total energy (kinetic energy 
plus potential energy) is positive. This criterion is different 
from Goodwin (1997), who marked all stars outside a given 



tidal radius as unbound. 

Fig. shows the evolution of the Lagrangian radii con- 
taining 10% and 50% of the current bound mass of the sys- 
tem and the virial ratio rj — — 2£'kin/£'pot of the bound 
particles. The constant mass radii and virial ratios before 
gas expulsion show that initially the system is indeed in 
virial equilibrium. When the superimposed gas potential de- 
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100 150 200 

dimensionless time t 

Figure 2. Time evolution of tlie mass radii (upper panel) and 
the virial ratio r] (lower panel) of the system shown in Fig. fll 



creases (t > 20), the mass radii increase rapidly and relax 
for t > toxp- This behaviour is also reflected in the virial ra- 
tio. The system achieves a more extended equilibrium state. 

The radial expansion factor of the cluster can be esti- 
mated as follows: In the adiabatic case (expulsion timescale 
is much longer than the c rossing time) R ■ M is constant 
and therefore ( |Hills 1980| ; |Mathieu 198^ the ratio of final 
to initial radius is 
Rf _ M^ _1 



with < e < 1. 



(2) 



On the other hand, if the expulsion time is short compared 
to the crossing time, we can apply conservation of kine tic en - 
ergy per parti cle du ring the ejection of the gas. Hills (198C) 
and Mathieu (1984) obtained 

1 = 27^ with i<e<l. (3) 

If e < 0.5 the final system is unbound. 

Fig. H shows the ratio of the final to the initial half- 
mass radii of the bound particles at the end of the simula- 
tions, compared to the theoretical predictions of Eq. Hand pi 
If stars are lost and the bound mass of the cluster is not 
conserved, Eq. H and H are not strictly valid any more and 
discrepancies to the analytic approximations occur. 

As expected, the simulations with long gas expulsion 
timescales fit well the solid curve, representing the adiabatic 
case. The faster the gas expulsion, the larger is the ratio of 
the final to the initial radius compared to the theoretical 
result. 

The models with fast gas expulsion follow the dashed 
curve well for high SFEs. For low SFEs, the ratio of final to 
initial radii is smaller than the analytical prediction, which 
emphasizes that the divergence for e = 0.5 does not oc- 
cur in numeric simulations: The final radius is decreased by 
excluding the unbound particles which are located preferen- 
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Figure 3. Ratio of the final to the initial half— mass radii ver- 
sus SFE of the simulations (model N3) for various expulsion 
timescales toxa 



tially at high radii. Additional, the outgoing particles reduce 
the total energy of the remaining system and leave behind 
a tighter bound core. 



2.3 Constraints on the Star Formation Efficiency 

After the gas expulsion the system of the remaining bodies 
relaxes again (Fig. 0). 

Fig. shows the ratio of the number of finally bound 
stars to the initial number of stars for various SFEs and gas 
expulsion timescales. The upper panel provides a resolution 
study of the runs Nl and N3 with 1000 and 4000 particles, 
respectively. Within the uncertainties they are indistinguish- 
able. Large dots show results from Lada et al. ( 1984 ) ob- 
tained from simulations with 50 (!) stars. As can be seen, 
the number of particles used does not influence the results. 

The lower panel compares the initially more concen- 
trated King model (N2) to the less concentrated one (Nl). 
We find that the curves of model N2 with f cxn > are shifted 



to lower star formation efficiencies or higher ratios of bound 
stars, respectively. This is due to the lower half-mass cross- 
ing time of the more concentrated cluster N2 (Table hi), con- 
firming that only the ratio of the expulsion timescale to the 
crossing time is important. Thus, more concentrated clusters 
have a larger chance to survive. In the case of instantaneous 
gas expulsion (icxp = 0) the models Nl and N2 yield the 
same curve. 

The ratio of bound to unbound stars gives a threshold 
for the SFE e necessary to form bound clusters. In case of 
instantaneous gas expulsion (see Fig. W) the curves are cen- 
tered around e = 0.45, which is somewhat less than the theo- 



retical limit 6 = .5 for bound clusters given by Hills (1980). 
Adams (EOOOl) recently gave analytic approximations for 
the dependency of the number of bound stars on the SFE in 
the case of instantaneous gas expulsion. For a SFE e = 0.5 
about 73% of the stars are kept, in good agreement to our 
results from Fig. (icxp = 0). However, our results show a 
stronger depe ndence of mass loss on the SFE e. Contrary to 
Adams (200C), star clusters with a SFE lower than e = 0.4 
are dissolved in our simulations. This discrepancy can be 
understood from the fact that Adams uses density distri- 
butions of gas and stars with very different concentrations. 
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Figure 4. Ratio of the number of bound stars to the initial num- 
ber of stars in the relaxed system after gas expulsion. Upper panel: 
models Nl and N3; lower panel: models Nl and N2; each symbol 
represents one run with given SFE and gas ex pulsion time; large 
dots are results taken from Lada et al. (L984), Fig. 2 therein 



Thus, even if the global SFE is small, the local SFE in the 
region of star formation could be as high as 90%, leading to 
a bound system. 

The number of finally bound stars increases with the 
gas expulsion time. In order for more than 50% of all the 
stars to remain bound the SFE must be equal to 45%, 30%, 
25% and 20% for gas expulsion times tcxp =0, 2, 4 and 10, 
respectively. The Galactic average SFE in giant molecular 
clouds is of the order of a few p ercent (Myers et al. 1986 
Williams fc McKee 1997|). Koo ([l999|) has observed SFEs 



up to 15% in the star forming-region W51B, maybe due to 
shock-interaction with a spiral density wave. Given these 
SFEs, unrealistic high gas expulsion timescales are required 
to obtain bound clusters. Only few clouds show SFEs up to 
30% or 40% (Lada 1992) and may stay bound. 



3 COMBINED N-BODY AND 
HYDRODYNAMICAL SIMULATIONS 

In Sect. B, the effect of gas removal is treated as a time 
variable external potential. To describe the physics more 
properly, we extend our simulations using smoothed parti- 
cle hydrodynamics (SPH, see Monaghan 1992 for a review). 
We are using an SPH-code with variable smoothing length, 
individual particle timesteps (Bate, Bonnell & Price 1995) 
and collisionless N-body particles, which was made available 
by Matthew Bate. 

The conversion from code units to physical units is the 



Table 2. Parameters of the initial configurations including gas 
dynamics (N-body & SPH) 



model 



rh 



N 



R 



SI 


0.7 


0.8 


2 X 4000 


0.76 


0.01 


1.0 


4.0 


S2 


0.7 


0.8 


2 X 11045 


0.76 


0.05 


1.0 


4.0 



rf^: half— mass radius; tc'. crossing time at half-mass radius; A^: 
number of particles; T: gas temperature; <5: numerical (Plummer) 
smoothing length (N-body part only); R: cu t-off radius; ^5: di- 
mensionless cut-off radius, see Bonnor ([19560; total mass (stars 
and gas) of all models was set to 1; all quantities are given in 
dimensionless code units (see text) 
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Figure 5. Evolution of the mass radii of a combined N-body and 
SPH simulation. At the end of the stability test the system is in 
equilibrium 



same as in Sect. 2.1. Additionally, we describe the internal 
energy of the gas with a dimensionless temperature that 
scales with T = G M fj,R~^ R^^'y~^ , where fi is the mean 
molecular weight, i?gas is the gas constant and 7 the constant 
adiabatic exponent. For an ideal gas with 7 = 3/2 and jj, = 
1.0, we have f = 3.5 10^ K. 



3.1 Initial Configuration and Models for Gas 
Expulsion 



Murray & Lin ( [1992D proposed as initial conditions pressure- 
confined protocluster clouds. As a reasonable initial config- 
uration we therefore adopt a pressure-c o nfined, isothe rmal 
Bonnor-Ebert (BE) sphere (Ebert 1955; Bonnor 1956), see 
Table |. 

To get a combined system of gas and stars, we add N- 
body particles with an equal density distribution, scaling 
the masses of the stars and the gas according to the SFE e. 
The velocity dispersion of the particles is chosen according 
to the temperature of the gas. 

Isothermal spheres extend to infinity. In order to get a 
finite configuration, the gas density is set to zero at an arbi- 
trary radius. To stabilize the gas sphere an external pressure 
is applied. This is not possible for particle systems. Therefore 
the velocity dispersion of the N-body particles is decreased 
and the system is allowed to relax until a stable configura- 
tion is obtained. 

Fig. U shows the evolution of one typical setup prior to 
gas removal. About 2% of the N-body particles are lost, but 
after some oscillations the main part achieves an equilibrium 
state. This configuration is then used as initial model for fol- 
lowing investigations. The different models used are shown 
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Figure 6. Time evolution of model S2 during the gas expulsion 
phase (dimensionless time 4). The SFE is e = 0.4 and central 
heating is applied. The plots show the N— body particles (small 
dots) projected onto the x—y— plane (only every 4th particle is 
shown); unbound particles are indicated by squares. The contour 
lines indicate the gas density in the x-y-plane in steps of 0.025. 



in Table H 

As the processes of gas expulsion are not well under- 
stood, we choose two simplified scenarios: In our first model 
we heat up the whole gas cloud. As a result, the gas starts to 
expand and is removed. Such a situation might be caused by 
stellar winds or ionizing radiation. In our second scenario we 
heat up a small inner core, creating an outward propagat- 
ing shock front disrupting the gas cloud. Such shock fronts 
may be generated by combined supernova explosions and 
winds from central high-mass stars. The formation of su- 
pershells and their ability to disrupt the cloud are discussed 
in var ious papers with regard to chemical self-e n richment of 



GCs ( [Morgan fc Lake 1989| ; [Brown et al. 199E|; [Parmentier 
et al. 11)99). Goodwin, Pearce fc Thomas (|2000) investigate 



single supernovae in gas clouds. In our simulations the heat- 
ing of the cloud, determining the gas expulsion timescale, 
must be sufficient to expel all the gas. 



3.2 Evolution of the Cluster During and After 
Gas Expulsion 

Fig. y shows the time evolution of a system with gas ejection 
in a supershell. Once the internal pressure of the expand- 
ing gas is smaller than the adopted external pressure, the 
outward propagating shell becomes unstable and develops 
substructures. At that stage, less than two crossing times 
after the gas removal started, the gas density in the cluster 
region is so low that its gravitational effect on the stellar 
component is negligible. We therefore remove the gas and 
follow the evolution of the stars alone. The globally heated 
cloud models show a very similar evolution and are treated 
in the same manner. 
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Figure 7. Ratio of the number of bound stars to the initial num- 
ber of stars in the relaxed system after gas expulsion; dashed lines 
indicate the results from Fig. W, model N3 



Fig. M corresponds to Fig. W in the pure N-body case. 
For comparison, the dashed curves reproduce the results of 
the pure N-body simulations. 

For instantaneous gas expulsion (all the gas particles 
are removed at once), corresponding to the case tcxp = 
in Sect. H, the simulations using the BE sphere as initial 
configuration are in very good agreement with the simula- 
tions using a King density distribution. The slightly higher 
crossing times of the BE sphere (see Table 0) may cause the 
small differences for small SFEs: the simulation was not run 
long enough for all particles to get unbound. We therefore 
conclude that the density distribution has no infiuence on 
the number of bound particles after gas expulsion, at least 
xp =0. 

The crosses in Fig. M, lower panel, show system SI where 
the temperature of the whole cloud was increased by a factor 
of 10 with respect to the equilibrium model. The number of 
bound stars increases slightly compared to the case with in- 
stantaneous gas expulsion. The diamonds in the upper panel 
show the results for a cloud centrally heated to T = 152 
(system S2). Compared to the first model SI, no significant 
differences are visible. Again, the number of bound stars in- 
creases slightly. However, the gas expulsion process in both 
cases is much faster than the timescales adopted in Sect. y. 
Therefore, also in more realistic cases, high SFEs are needed 
to sustain a bound star cluster. 

One way out may be a collapse of the star cluster before 
the gas is complet ely ex pelled, leading to a high er "effective 
SFE" . Lada et al. ( |l984| ) and Verschueren ( |l990| ) proposed a 
low or zero initial velocity dispersion to explain the coUaps. 
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Saiyadpour, Deiss & Kegel (1997) considered the effect of 



dynamical friction on the stellar cluster. 

We implement the first approach by setting the initial 
velocity dispersion of the stars to zero and adopting gas ex- 
pulsion by heating the whole gas cloud. The cluster virializes 
within a small radius and only few stars that gain velocities 
higher than the escape velocity of the cluster are ejected. 
For a wide SFE range the percentage of bound stars at the 
end of the simulations is nearly constant and higher than 
80% (Fig. |7|, upper panel, triangles) . Even our run with the 
lowest SFE e = 0.02 leads to a remaining bound system. 
We conclude that this scenario can easily explain the for- 
mation of bound clusters. However, if the gas expulsion is 
delayed one dynamical time scale or longer, the stars virial- 
ize in a smaller volume. We then basically get back to the 
case of an initially virialized system with a higher local SFE 
than in the beginning. The system may break up again and 
the number of bound stars will be less than in the scenario 
with low velocity dispersion, but higher than in the initially 
virialized case. 
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4 CONCLUSIONS & OUTLOOK 

N-body and combined N-body & SPH calculations to in- 
vestigate the influence of the residual gas expulsion on the 
stellar part in star forming regions are presented. 

We show that in the case of instantaneous gas expulsion, 
clusters with SFEs greater than e = 0.45 can keep more 
than 50% of the initial stars. Clusters with SFEs less than 
e = 0.40 are dissolved. Different concentrations of the initial 
models show no effect at all on the number of bound stars 
if the gas is expelled instantaneous. 

We confirm that gas expulsion timescales which are sev- 
eral times longer than the crossing time of the GC can de- 
crease the SFE needed to sustain bound clusters consider- 
ably. 

However, our simulations including the proper dynam- 
ics of the residual gas show that in order to destroy the whole 
cloud by global heating or by supershells the gas expulsion 
must take place on a short timescale, requiring a high SFE. 
Only few star forming regions show such high SFEs. 

We demonstrate that models with stars having an initial 
zero velocity dispersion lead to a compaction of the cluster 
and can explain bound systems even in low SFE regions: For 
SFEs as low as e = 0.15 more than 80% of the stars stay 
bound. Bound systems are obtained even with SFEs lower 
than 10%. For future investigations it is essential to know 
the velocity dispersion of newly born stars in clusters. 
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